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Abstract 

We propose a symmetrized form of the softened gravitational potential which is a natural extension of the Plummer 

potential. The gravitational potential at the position of particle ; (jc,, y,, z,), induced by particle j at (xj, jj, Zj), is given 

by: 

Gm i 



h/ + e;2 + e/|i/2' 



where G is the gravitational constant, nij is the mass of particle j, rtj - |(jc,- - xy)^ + (j; - yj)'^ + (z; - Zj)^\^^^ and e; 
and e, are the gravitational softening lengths of particles / and j, respectively. This form is formally an extension of 
the Newtonian potential to five dimensions. The derivative of this equation in the x,y, and z directions correspond to 
the gravitational accelerations in these directions and these are always symmetric between two particles. 

When one applies this potential to a group of particles with different softening lengths, as is the case with a tree 
code, an averaged gravitational softening length for the group can be used. We find that the most suitable averaged 
softening length for a group of particles is (ef) = 2^ mje//M, where M = YIj rrij and N are the mass and number of 

all particles in the group, respectively. The leading error related to the softening length is O [Yj rjd(e/)/ri/), where 
rj is the distance between particle j and the center of mass of the group and die/') - e/ — (e/). Using this averaged 
gravitational softening length with the tree method, one can use a single tree to evaluate the gravitational forces 
for a system of particles with a wide variety of gravitational softening lengths. Consequently, this will reduce the 
calculation cost of the gravitational force for such a system with different softenings without the need for complicated 
forms of softening. We present the result of simple numerical tests. We found that our modification of the Plummer 
potential works well. 
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1. Introduction 



(iPlumme 



I 



AarsethI (Il963h introduced a Plummer potential 
19111) for gravitational A^-body simulations 



in order to eliminate the divergence of the potential that 
occurs when two particles undergo a close encounter. 
With a Plummer softening, the modified gravitational 
potential between two particles / and j is given by 

Gmj 



mass of particle j and the gravitational softening 
length, respectively. The increase in the calcula- 
tion cost due to this modification of the potential is 
quite small. Thus, it is widely used in simulations 
of self-gravitating systems and the GRAPE (GRAv- 
ity PipE) series also adopts this form of the gravita- 
tional potential dSugimoto et all Il99rt 'ito et a l.L Il991 
Okumura et al.l[l993; Makino et al.]l997; Kawai et al 



\rii 



(1) 



2OOOI: iMakino et all [20031: iKawai & Fukushigel 12006 : 



Makinoi 120051 l2007h . Due to the wide adoption of 



where G, 
stant, the 



rij, nij, and e are the gravitational con- 
distance between particles / and j, the 
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the Plummer potential, there have been several stud- 
ies on the c hoice o f the optimal gravitational s often- 
ing length d Merritd 1 19961: lAthanassoula et all I2OO0I: 
iDehnen,200lT 
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In simulations of self-gravitating systems, it is nec- 
essary to deal with a wide range in both space and 
time. For example, one needs to cover the range of 
~ 100 Mpc (for the tidal field) to ~ 10 pc or less (for the 
structure of the interstellar medium) in galaxy forma- 
tion simulations . Usually, the outer region of the target 
galaxy is represented by massiv e particles in order to re- 
duce the calculation cost ('e.g.. lNavarro & Bena . Il99ll : 
Katz & WhiteLll993l) . In simulations of the coevolution 
of the galactic center and star cluster, it is necessary to 
use different masses for different components, since re- 
quired accuracies and s patial resolutions are diff'erent 
(IFuiii et all 120071 BOIOl) . Thus we need to use differ- 
ent gravitational softening lengths for different types of 
particles. 

If one uses the direct summation method to calcu- 
late the gravitational force, it is possible to use an ar- 
bitrary symmetric form for the softening length; for ex- 
ample, [(e,-2 + e/)/2]'/2, (e,. + e^.)/2, or max(e„ e,). On 



the other hand, when one uses the tree method (lAppelL 
19851 : Barnes & Hutl 1 19861) , one needs to use separate 
trees for each of the different groups of p articles with the 
same gravitational softening length (e.g. lSpringel et al. , 
20011) . since otherwise there will be an error in the force 



calculation of the order e . If one were to calculate 
the force from the center of mass of two particles with 
widely different gravitational softening values, it is not 
clear how one would take into account the gravitational 
softening length. However, if one uses different trees for 
particles with different gravitational softening lengths, 
th e calculation cost i n crease s by a large factor. 

Hernquist & Kata ( Il989b proposed a different ap- 
proach for softening. They assumed that each parti- 
cle has a finite size and mass distribution given by a 
function with a compact support, e.g. a spline func- 
tion. The advantage of this gravity kernel is that it be- 
comes a purely Keplarian potential for distances larger 
than the particle size. Thus, for particles separated by 
distances larger than their sizes, one can use the center- 
of-mass approximation or the multipole expansion of 
a purely Keplarian potential. However, spline func- 
tions require a few conditional branches and the eval- 
uation of a polynomial. These operations lead to a non- 
negligible slowdown in the force and potential calcula- 
tions on modern CPUs/GPUs which have SIMD execu- 
tion units. As a result, the calculation cost of a spline 
function is higher than that of the simple Plummer soft- 
ening. Moreover, the tree approximation cannot be used 
for particles within the softening length. Thus, if there 
are many particles within the softening length, they can 
cause a large increase in the calculation cost. 

In this paper, we propose a symmetrized gravitational 



potential between two particles with different soften- 
ings, which is a natural extension of the Plummer po- 
tential. We derive the multipole expansion of a group 
of such particles with this symmetrized potential. We 
then show that, in terms of accuracy, it is the best to use 
the second moment of the softening lengths of the par- 
ticles in order to obtain the averaged softening length 
for the group. The multipole moment up to the second 
order has a simple form, which is easy and inexpensive 
to calculate. 

This paper is structured as follows: in section |2] we 
present a symmetrized gravitational potential for a pair 
of particles with different gravitational softening lengths 
and then introduce an averaged gravitational softening 
length for a group of such particles. The use of this 
potential with the tree method is also discussed in the 
section. Results of numerical tests are given in section 
|3] A summary and discussion are shown in section|4] 

2. Symmetrized Plummer Potential and its Applica- 
tions to the Tree Code 

2.1. Symmetrized Pairwise Potential 

We modify the Plummer potential (Eq. [TJ so that it is 
symmetric between two particles even when their grav- 
itational softening lengths are different. Here, we con- 
sider the following form: 



Gmj 



+ 6,.2 + 6/|l/2- 



(2) 



Formally, this form can be regarded as the 1 /r potential 
in the five dimensional space (x,y, z, e,, e,). This form 
has been used before to model the gravitational inter- 
action between_gak}a;_particks_wiA different size and 
mass JWhiteL 1 19761: lAarseth & FallL Il980h . When we 
take the gradients of Eq. |2] in the x, y and z directions, 
we obtain the gravitational accelerations, i.e. 



^IJ,X 



dx 

Gmjixi - Xj) 



(3) 



This gravitational acceleration satisfies Newton's third 
law. 

2.2. Multipole Expansion of the Potential of a Group of 
Particles with the Symmetrized Plummer Potential 

In this section, we derive the multipole expansion of 
the gravitational potential for a group of particles with 
the symmetrized Plummer potential by introducing an 



averaged softening for the group of particles. Using Eq. 
|2] we obtain the most suitable form of the averaged soft- 
ening which minimizes the potential/fo rce error. This 
form can easily be used with a tree code (iBarnes & Huj 
19861) . 



We consider the gravitational potential at r, due to a 
group of particles j - I, . . .,N with the total mass M. 
The three-dimensional position of particle j is expressed 
as rj. Here we choose the center of mass of particles j 
as the coordinate center. By taking the summation over 
j, the potential at r,- is 



<Pi 



■Z 



Gnij 



|r,-/ + e,.2 + e/|i/2' 



(4) 



where r,j = r, - rj. 

We consider the following form as the monopole po- 
tential of the group: 



GM 



|r,.2 + e,2 + £/|i/2' 



(5) 



where fi/ is an arbitrary form of the averaged softening 
length for the group of particles j. As shown below, 
in terms of accuracy, we can obtain the most suitable 
form of the averaged softening length for the group of 
particles. 

By using Sf, we can rewrite Eq. |4]as follows: 



■Z 



Gmj 



|(^._^.)2 + e,2+ £.2+^(^.2)1 1/2 



(6) 



where (i(e/^) = e/ - £/. The Taylor expansion of Eq. 
|6]up to the second order of rj and d{ej^) is 
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(7) 



7; = (r,2 + 6,2 +£/)'". (8) 

The zeroth term of the Taylor expansion is 0', itself (Eq. 
|5]l. The term of the order of (r, ■ ry) vanishes by defini- 
tion. Since d(ej^) - e/ - fi/^, the first order term of 
d{ej-) vanishes if we adopt the second moment of e, as 
the averaged softening length: 



s/ 






J "J 



(9) 



Hence, choosing the second moment of Cj to be the av- 
eraged softening length is the most favourable option; 
using other forms such as (e,- -i- e/)/2, max(e,-, e,) and 
min(ei, e,) result in a second-order error term. 

For the convergence of the multipole expansion, it 
is necessary to take into account not only the ordinary 
three-dimensional distance but also the distribution of 
the softenings, £j. We will discuss this point in the next 
subsection. 

2.3. The opening criterion 

In the tree method, the force and potential of a group 
of particles are approximated by that of t he correspond- 
ing tr ee node. The opening criterion ( IBarnes & Hut , 

' 19861) . 

w 
e>-;, (10) 

d 

is widely used to determine whether the multipole ap- 
proximation can be used or not. Here, 6 is called the 
opening angle and controls the force/potential accura- 
cies, w is the size of the tree node and d is the distance 
between the barycenter of the tree node and the posi- 
tion where we want to calculate the force and poten- 
tial. Since the multipole expansion of the 1 /r poten- 
tial around r + dr is a series of dr/r where dr ~ w 
and r ~ d, Eq. \W\ corresponds to the convergence 
condition for the multipole expansion of the 1/r poten- 
tial. For the Plummer potential, the Plummer distance, 
{r- + e,2)'^2, is adopted instead of r. Typically, a value of 
= 0.5 ~ 0.75 is used. Note that, unless 6 < 1/ V3, the 
Barnes & Hut's opening criterion causes unbound errors 
in force calculation, resu lting in "detonating galaxies" 
JSalmon & Warrenill994l) . They studied other opening 
criteria so that one can avoid these large errors. Here, 
we use Eq. [10] with sufficiently small 0. Hence our re- 
sults, which will appear in section [331 are free from this 
problem. 

As shown in section 12.21 the potential induced by a 
group of particles with the symmetrized form of the 
potential is expanded by two parameters, rj and c/e/. 
Thus, the following two constraints together should sat- 
isfy the conditions for convergence: 



1 > — > 
' R 



and 



d(ej^) 
R2 ■ 



(11) 



(12) 



Here, r/ and rj' are the tolerance parameters for the mul- 
tipole expansion of the symmetrized Plummer potential. 
By replacing rj with w and dej^ with Cnax^—emm^, where 



fmax and e^in are the maximum and minimum softening 
lengths in a tree node, respectively, we obtain 



77> 



R' 



and 
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(13) 



(14) 



for the convergence conditions of the multipole expan- 
sion. 

Note that Smax^ - fmin^ corrcsponds to the size of 
the "box" in the fifth dimension. In our current im- 
plementation, the tree structure is still an oct-tree in 
three dimensions. Thus, if a node contains two parti- 
cles which have widely different gravitational soften- 
ing lengths, its "size" can be dominated by the soften- 
ing if its physical size is small. When we calculate the 
force from this node to a nearby particle with a distance 
^ [(fmax^ - fmin^)/??']'^^. this nodc is always opened up 
even if the physical size of the node in three dimensions 
is much smaller than the distance. This behavior could 
cause an unnecessary increase in the calculation cost. In 
principle, we can easily solve this problem by making a 
tree structure in four-dimensional space, but so far we 
have not used such treatment. 

In the tree-with-GRAPE method (M akinoi Il99lb . a 
similar opening criterion, a set of Eqs [13] and [141 can be 
used with two modifications based on an argument for 
the convergence condition of the Taylor expansion. The 
modifications are as follows; First, we employ the min- 
imum node-to-node distance in the ordinary three di- 
mensions as the separation distance, r,, in Eq. [8] Then, 
we use the minimum gravitational softening length in 
the particles which share an interaction list, instead of e, 
in Eq. [8] These modifications give a sufficient accuracy 
when used with the tree-with-GRAPE method. 

2.4. Calculation of averaged softening 

The averaged gravitational softening length for each 
tree node can be calculated in the tree building phase 
with a recursive operation that is the same as that for the 
multipole expansion. The averaged gravitational soften- 
ing length for a parent node, (e^), is 



<fp'> = 






(15) 



where Nc is the number of child nodes, M^x and {ec,k^} 
are the mass and the averaged gravitational softening 
length of the k-th child node, respectively, and Mp is the 
mass of the parent node. 



2.5. Higher order expansion and application to FMM 

In this section, we discuss the application of the sym- 
metric softening to the fast rnultipol e method (FMM) 
dGreengard & RokhlinL Il987l: IZhaol Il987t JDehneni 
2000). UnUke the tree method, in FMM, multipole ex- 
pansions are first translated to local expansion, which 
is then evaluated at the positions of multiple particles. 
Thus, the calculation order in FMM is proportional to 

The calculation cost in FMM depends on the 
expansion order of the multipole moments in the 
local/external particle s. In the standard FMM 



dGreengard & Rokhlinl 119870 . spherical harmonics is 
used to express expansion and thus the calculation 
cost is Oip'^N), where p is the order of the multi- 
pole expansion. JElliott & Board! ( 1 19961) successfully 
reduced the calculation cost in FMM from 0{p'^N) 
to 0{p{logpy^~N) by using a fast fourier tra nsforma- 
tion a nd optimal choice of the tree level. JDehnenl 
(I2OO0I) proposed a Cartesian FMM, however, its cal- 
culation cost is estimated as 0{p^N). Again this 
cost can be reduced to 0{p^^^{log p^^^N) by following 
Elliott &Boardl(ll996l) . 

With the 1/r potential in the five dimensions pre- 
sented here, one cannot use a spherical harmonic func- 
tion for the expansion. Thus, it is necessary to use a 
straightforward expansion such as what is used in the 
Cartesian FMM ( lDehnenil2000l) to implement our sym- 
metrized Plummer potential. In the multipole expan- 
sions, because it is necessary to expand the fourth or 
fifth dimension for the local and external field, the cal- 
culation cost is O(p^N) and the reduced calculation cost 
is Oip^ilogp)^^^ N). When the multipole expansion is 
applied, the softening lengths of the particles are usually 
smaller than the size of the node and therefore one usu- 
ally needs to take into account only the first cross term in 
Eq. [7] In this case, the calculation cost is similar to that 
of the Cartesian FMM. Moreover, gravitational soften- 
ing is used in simulations where the required force accu- 
racy is not very high. Therefore the multipole expansion 
can be truncated in relatively low orders. Thus, the in- 
crease in the calculation cost is not so large if one incor- 
porates the symmetrized Plummer potential into FMM. 



3. Tests 

3.1. Error in Acceleration 

Here, we present the measured acceleration error of 
the force calculated using a tree and our symmetrized 
form of the Plummer potential. We take two groups of 



5 X lO'* particles, each with a different mass and soft- 
ening length. The particles are randomly distributed 
within a uniform sphere of radius L = 1. Table [1] shows 
the different combinations of masses and gravitational 
softening lengths used for this test. The accuracy of 
the accelera tions depends s trongl y on the distribution 
of particles dBarnes & Hull Il989b . However, for sim- 
plicity, we consider only a uniform sphere here. In this 
test, we adopt the center-of-mass approximation for the 
force calculation and we use the set of Eqs[T3]andfT4lfor 
the opening criterion. For this criterion, the tolerance is 
always kept the same and is represented by 9 instead of 
7/ and T]' for simplicity. 

We measure the mean acceleration error as a function 
of the opening angle for the three different mass ratios: 
1:1, 1:8 and 1:64 (see Table [TJ. The mean acceleration 
error for each case is given by the following equation: 



1 /. \ai{e) 



■ «i.O 



(16) 
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where a,(0) is the acceleration vector of particle / es- 
timated with tolerance 6 and a,_o is that evaluated for 
9 = 0. We also measure the mean acceleration error 
which excludes the error due to d{e^). To do this, we 
build two trees, where each tree consists of particles 
with the same mass and softening. We then calculate 
the accelerations of the particles by adding the forces 
from the two trees so that we obtain the mean accelera- 
tion error without the contribution of d{£^). 

The mean acceleration error as a function of 6, for 
each of these methods, is shown in Figure [1] We found 
that there is no clear difference in the mean accelera- 
tion errors between different methods. We also found 
that the mean acceleration error is approximately pro- 
portional to 6^ . This is in good agreement with previous 
results obtaine d by the tree code with the usual Plum- 
mer potential JHernquistl 1 19871: iBarnes & Hul Il989 : 
Makinoi Il990l) . This means that our adopted opening 
criteria work correctly. 

3.2. Calculation Cost 

To evaluate the calculation cost, we measured the av- 
erage number of particle-particle interactions and the 
average number of total interactions. Figure |2] shows 
the these numbers as a function of 6. We used the same 
distribution of particles as we used in section 13.11 The 
dotted curves are the calculation cost for the case of 
two separate trees for particles with different softenings. 
The calculation cost of our method is about one half of 
that of the calculation with two separate trees. The av- 
erage number of the interactions is roughly proportional 
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Figure 1 : Mean acceleration error as a function of the tolerance pa- 
rameter 0. Curves with crosses, circles and triangles are the results 
for mil mi = 1,8 and 64, respectively. The solid and dashed curves 
are the results obtained by the tree which uses our averaged softening 
length ("This method"), whilst the dashed curves are those obtained 
from using two separate trees, where each tree consists of each parti- 
cle group ("Separate"). 



to ^ for both methods. This dependence is again close 
to those in previous works {Makino^ il99Q\ 

3.3. Energy Error in Time Integration 

Here, we present the result of test calculations of a 
cold collapse. We create a uniform sphere with the 
radius L and mass M using a random distribution of 
1024 particles, where the system of units we use is 
L = M = G = 1 . The sphere consists of two groups 
of 512 particles. The mass and softening length of par- 
ticles in these two groups are varied and the values used 
are given in Table |2] The initial virial ratio of the sys- 
tem was set to 0. 1 and the initial particle velocities were 
drawn from a Gaussian distribution. 

The numerical simulations were performed using 
ASURA (Saitoh, in preparation). ASURA adopts the tree- 
with-GRAPE method to calculate the gravitational force 
and potential ( Makino, 1991) . For this test, ASURA only 
uses the center-of-mass of the tree nodes to calculate 
the force, with their averaged gravitational softening 
lengths determined according to Eq. |9]and the tolerance 
is set to 0.5. The maximum nu mber of partic l es which 
shares the same interaction list dBarnesl 1 1 990l : iMakino , 
Il99lh was set to n„ = 8 and n^ = 128. Time-integration 
scheme used is an ordinary leap-frog method. 



Table 1 : The summary of the masses and gravitational softening lengths of the two types of particles used in our test simulations, m i and ei refer 
to the mass and gravitational softening length for the less massive particles, whereas m? and £2 are those of the massive particles. The number of 
particles in each group is 5 X 10"*. 



Mass ratio 


nil 


ei 


1112 


£2 


1 : 1 

1 :8 
1 : 64 


l.OOx 10-5 
2.22 X lO-*- 
3.08 X 10-^ 


6.79 X 10-' 
4.11 X 10-^^ 
2.13x10-^^ 


l.OOx 10-5 
1.78x10-5 
1.97 X 10-5 


6.79 X 10-' 
8.22x10-^^ 
8.51 X 10-^^ 


Table 2: 


The same as table [T]but for the test calculations of cold 


collapse. 


Mass ratio 


OTi 


ei 


»Z2 


£2 


1 : 1 
1 :8 
1 :64 


9.77 X lO-'' 
2.17x10-'* 
3.00x10-5 
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1.89 X 10-2 
9.79 X 10-^^ 
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1.92x10-^^ 
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Figure 2: Number of the average interactions as a function of 8 with 
L = 1 and mi/oii = 64. Curves with symbols are the number of the 
mean particle-particle interactions, whereas curves without symbols 
are the sum of the number of the mean particle-particle and mean 
particle-node interactions. 



Figure[3]shows the relative error in the total energy of 
the system, AE, at time t - I (which is approximately 
the time of maximum collapse) as a function of the size 
of the time-step. Here AE is defined as 



AE 



£(f = 0)-£(f = 1) 



Eit = 0) 



(17) 



where E{t) is the total energy of the system at t. 

From figure[3] we can see that the mass ratio between 
groups of particles does not affect the relative error. For 
small time-steps (< 2-^), the relative error of the tree 
method is larger than that of the direct method. This 
is because the error in the calculated force from using 
a tree dominates the integration error The contribution 
of the tree approximation to the relative error is smaller 
for large rig, since a larger numbe r of nearby par ticles is 
included in the direct calculation ( lBarneslll990l) . 



4. Summary and Discussion 

In this paper, we introduced a natural symmetrization 
of the Plummer potential. This symmetrization is quite 
simple and hence, easy for anyone to implement into 
their own code, providing their code adopts the tradi- 
tional Plummer potential for gravity. 

By applying this potential to a group of particles, we 
derived an averaged gravitational softening length for 
the group. Consequently, we can use this potential in 
the tree method. The multipole expansion of the sym- 
metrized Plummer potential allows us to use a single 
tree to calculate forces and potentials for a system of 
particles that have different softening lengths. Thus the 
calculation cost of our method is less than that of previ- 
ous implementations. 
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Figure 3: Dependence of the relative energy error (at time f = 1) 
on time-step size. Crosses, circles and triangles are the results for 
m2/fni = 1,8 and 64, respectively. Solid, dotted, and dashed curves 
are the results for rig = 128, fig = 8, and direct calculation, respec- 
tively. 



Since the modification is quite simple, the modified 
Plummer potential can be easily implemented in mod- 
ern GRAPEs (i.e., GRAPE-7[]and GRAPE-DrEI) as 
they can modify the force law via software libraries. 
However, for older GRAPEs, the gravity pipelines are 
implemented on the hardware layer and so it is impossi- 
ble to use this modified potential. The latest version of 
the GRAPE-DR library supports this new symmetrized 
Plummer potential. 

Pure software libraries, i.e. those that run on typ- 
ical CPUs without special hardware, can easily in- 
clude our calculation model. Phantom-GRAPE (Ni- 
tadori et al, in preparation) is an optimized library 
that utiUzes the SIMD instructions of X86 CPUs and 
has application inte rfaces comparable with GRAPE-5 
dKawai et al.L 120001) . For this software, we can eas- 
ily modify the gravity kernel in the library. We used 
this modified Phantom-GRAPE for the tests described 
in Section[3] 

The implemen tation of our modified Plummer poten- 
tial on GPGPUs ( lNakasatoll2009HHamada et allboogl) 
is also straightforward. 



GRAPE-7 employed field programmable g ate array (FPGA) and 
constructs the gravity pipelines on runtime JKawai & Fukushigg . 
I2006h . 

^GRAPE-DR uses the original programmable SIMD processors 
iMaldnoLl20o3 : lMakino et alir2007h . 
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